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Abstract 

At the 1999 Game Developers Conference, Thomas Jakobsen presented a semi-implicit rigid 
body physics simulation. The simulation could integrate stiff systems at the time steps used in 
games without the need to compute expensive Jacobian matrices and had very simple 
constraint implementations. However, the simulation was limited because it did not include a 
velocity state and it did not integrate the rotational equations of motion, necessitating large 
numbers of particles and distance constraints to simulate a single rigid body. It was also 
difficult to set up pin constraints, common in rag-doll simulations. 

By extending the simulation to six degrees of freedom through the Verlet-like integration of a 
state quaternion, joint constraints became simple and straightforward, though a slightly more 
complicated constraint equation was required. Two methods for implementing the additional 
degrees of freedom are presented; one that uses the current quaternion and last quaternion as 
the system states and one that uses the quaternion and the angular rate. The advantages and 
disadvantages of both methods are discussed. The equations used in pin and angle 
constraints and the method for resolving collisions are presented. This simulation was 
developed for the game MX Unleashed, published by THQ. 


1 Introduction 

The decision to use the simulation method described in this paper was made by discovering 
the disadvantages of most of the simulation techniques that have already been published. 
Initially, a simulation using explicitly integrated equations of motion was developed. This 
simulation had the familiar stability problems when trying to integrate stiff systems. Constraints 
were formulated using the methods described by [Witkin]; write the constraint equation, 
differentiate it twice and solve for the constraint forces. The constraint solver worked, but some 
of the constraints had mathematically indeterminate solutions at specific orientations of the 
rigid bodies that had to be carefully avoided. Collision response was handled through the use 



of impulses [Baraff] since penalty methods could not be made stiff enough. Unfortunately, the 
rigid body objects would pop once they had come to rest on a surface. 

Implicit integrators were used in order to handle stiffer systems. However, it was next to 
impossible to derive an expression for the Jacobian, or force derivative matrix, for some of the 
most computationally intensive, non-linear force models, like the tire model in a car simulation. 
Instead, the Jacobian was computed numerically by calling these force models a couple of 
times with small perturbations in each of the states. The cost of doing this, combined with the 
cost of solving a large set of linear equations, became prohibitive, especially when trying to 
simulate many four wheeled vehicles. 

The simulation presented by [Jakobsen] was implemented next. The semi-implicit Verlet-like 
integrator permitted the use of springs stiff enough to handle resting contacts without the need 
to compute a Jacobian and the simple constraints always worked. However, since the 
simulation did not include a velocity state and did not integrate the rotational degrees of 
freedom, it had to be extended to be practical in our application. The final result met all of our 
needs and is described herein. 


2 The Particle 

The basic particle is represented as a point of mass, m, that has a position, x(0 and a velocity, 
v(0, that are functions of time. 


3 The Rigid Body 

For the purposes of this simulation, a rigid body is approximated by a particle with a mass 
distribution, or inertia, I, that has orientation, represented by a quaternion, q(r). 
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(3.1) 


I is not a function of time. It is computed once in the rigid body coordinate frame. This saves 
the cost of rotating the inertia matrix into the world coordinate frame each time step. The 
complete inertia matrix, I, can be used in the equations of motion that follow. However, by 
setting all of the products of inertia, or off-diagonal terms, to zero, the matrix can be replaced 
with a vector. Equation (3.2), which reduces the number of calculations that must be performed 
throughout the simulation. 
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4 Equations of Motion 


To animate the particle, the simulation presented by [Jakobsen] was used. However, since the 
simulation only integrated the positions of particles, sets of particles had to be connected 
together with distance constraints in order to simulate a rigid body. Consequently, it was more 
expensive than including the rotational degrees of freedom. 

The [Jakobsen] integrator was modified slightly to introduce a damping term: 

x(r- Ar) = x(r) 

The i^veiocity term can be used to introduce artificial damping into these equations of motion. 
Normally, ^vehaty equals one, and this equation is identical to [Jakobsen]. When the damping 
constant is zero, the equations of motion are reduced to F = mv; as soon as the external force 
is removed, the body stops moving. 

This method was extended to integrate the rotational degrees of freedom: 

6q = Slerp(q,q (r)q(r-Ar),^ 

angular rate ^ 

q'{t + At) = 6qq{t) 

q(r + Ar) = q(r) + q(q'(r + Ar),i 'T(r))Ar^ 
q(r-Ar) = q(r) 


where 
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Like ^velocity in the translational equations of motion, ^anguiarrate normally equals one. If the angular 
rate damping equals zero, the body will stop rotating when the external torque, t(0, is 
removed. \i ianguiar mte is always one, the Slerp can be skipped and 6q = q(r) q(r-Ar). The Slerp, or 
spherical linear interpolation function, along with some of the other quaternion math operations 
are listed in Appendix A. Equation (4.3) is reproduced from [Refson]. 



Initially, this worked well. However, by using the current state and the last state in the 
equations of motion, any force that is a function of the velocity or angular rate, like a damping 
term, must use a finite difference method to compute them. 
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These calculations can be very noisy and the problem is exacerbated when constraints come 
into play and instantaneously change the particle states. (Please note that e> = q{t) (Hworid q(0 
simply rotates the angular rate from world to body coordinates.) 

This problem is resolved in the translational degrees of freedom by replacing the x{t-At) state 
with velocity, \{t). The new equations of motion become: 
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[Garberoglio] provides a description of a Verlet-like integrator for quaternions. Like the 
translational equations of motion, the q(r-A0 state is replaced with body angular rate, e>(0- 


q(r + Ar) = q(r) + q(q(r),e>(r))Ar + ^q(q(r),e>(r),l^'T(r))Ar^ 
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Although this second set of equations provides smooth velocities and angular rates, the 
constraints must now deal with these new states. In the first set, referred to as the velocity-less 
equations of motion, although constraints can set the velocity and angular rate by setting 
x{t-At) and q(r-A0, they do not have to. At a minimum, they can simply reposition or reorient the 
body. Therefore, the set of equations used depends on the application. If the application does 
not generate forces that are a function of the rates, the velocity-less equations of motion can 



be used to take advantage of the very simple and fast constraint equations. Equation (4.8) is 
derived from [Refson]. 


5 Collision Response 

The simulation uses the [Baraff] method for resolving collisions between two rigid bodies. 
However, instead of computing the contact forces for resting contacts, which can be expensive 
when there are multiple contact points, an easier and faster penalty method is used. This 
penalty method can be problematic when using explicit integration, but the semi-implicit 
integrator permits the use of stiff springs. 



Figure 1 - Collision Geometry Output from Collision Detection Algorithm 


As illustrated in Figure 1, the collision detection algorithm returns a set of contact points and 
surface normals. For every contact point, c„ and surface normal, ft,, the following set of 
calculations is performed. First the velocity at the point of contact is computed. 

r,“=c,-x“(<) 

<„,.=q"(')“"(')r(') (5.1) 

v“ , ,= xr." 

contact \ / world i 


If a second rigid body is involved in the collision, the second body’s contact point velocity is 
computed using Equation (5.2), otherwise it is set to zero. 



(5.2) 


The relative velocity is simply: 

v" = v“ - V* 

relative contact contact 

V* = — 

relative relative 


(5.3) 


If the relative velocity along the surface normal is small, the contact point is considered a 
resting contact. 
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In which case, a penalty force is calculated so the two objects can rest on one another. 
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e, is the point on the surface of the rigid body that also lies on the ray between the center of 
gravity and the contact point, c, (Figure 1). If a second rigid body is involved in the collision, this 
force is nof applied to the second body. Instead, a second penetration force is calculated using 
a different set of spring and damping constants and both forces are scaled in relation to their 
masses, as presented in Equation (5.7). These forces are then applied to their respective rigid 
bodies. 
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A potential problem with this approach to collision resolution is that the resulting motion of the 
rigid body or bodies varies with the number of contact points returned from the collision 
detection algorithm. More points generate more resting and friction forces so ksprmg must be 
lowered to compensate. The collision detection algorithm used with this simulation always 
returns the same number of contact points for any particular collision orientation so this was 
not a problem. 

If the penetration speed, Sreiative, is large, impulses are used to instantaneously change the 
velocity and angular rate of the objects involved in the collision. Other than minor changes in 
notation. Equation (5.9) is presented unchanged from [Barraf] page D47, Equation (8-18). 
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s is the coefficient of restitution. This impulse is then used to change the velocities and angular 
rates of the two colliding bodies. 
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One may question the need for Equation (5.8) if the inertia matrix, r\ were transformed into 
the world coordinate frame for each of the rigid bodies. This savings comes at the cost of 
transforming the inertia matrix, I, into the world coordinate frame and then finding its inverse, 
which is more expensive than four vector rotations. Some simulations transform the inertia 
matrix into the world frame after the equations of motion are integrated. The problem with this 
approach is that I and its inverse, r\ would have to be re-transformed every time the 
orientation of the rigid body is updated by a constraint. 


The final step of the collision response is to compute a friction force using a simple kinetic 
friction model. This model transitions into a viscous friction model at low sliding speeds and as 
a result, ignores the effects of static friction. This was suitable for our purposes, but any friction 



model can be used. The friction coefficient, iikinetic, changes based on which surface the rigid 
body is sliding on. 
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6 The Distance Constraint 
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Figure 2 - Particle Distance Constraint 


The distance between two particles can be kept relatively constant through the use of the 
distance constraint (Figure 2). This is very useful in rope, cloth, water surface and soft body 
simulations. The simple distance constraint presented by [Jakobsen], page 388, assumes that 
the masses of the two particles are equal, which is generally true. However, to include the 
effects of mass, it is simply a matter of repositioning the particles as a function of their mass 
ratios. 
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If the masses are identical, the mass ratio terms can simply be replaced with !4. A constraint 
strength term, kstrength, which can range from zero to one, is introduced to soften the distance 
constraint. This allows the rope, cloth or soft body to stretch more than is possible by only 
lowering the number of constraint iterations. 

In games, the results are more than satisfactory. However, by ignoring the velocity states, 
residual motion is left in the particles that may not be desirable. Differentiating the constraint 
equation. Equation (6.1), yields Equation (6.3); the velocity difference between the two 
particles cannot move the particles closer or farther apart. 

c = (v‘'(r)-v*(r)^-6x = 0 (6.3) 
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If needed, these corrections can be made in the velocity-less equations of motion by replacing 
the expression for 6v in Equation (6.4) with the one in Equation (6.5) and then setting the 
previous positions of the two particles using Equation (6.6). However, if the velocity-less 
equations of motion are suitable for the application, these corrections are most likely 
unnecessary. 
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7 The Pin Constraint 

The distance constraint becomes a pin constraint when it is applied to rigid bodies. It connects 
two bodies at arbitrary points in or on the surface as illustrated in Figure 3. 



The vector between the two pin points on the two bodies is first calculated: 

( 7 - 1 ) 


The problem is how to divide 6x between the orientations and positions of the two bodies and 
in what proportions. Fortunately, the impulse equation from [Baraff], Equation (5.9), was re- 
applied to arrive at the solution: 
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If the velocity states are being used, the quaternions are not updated at this point because 
they are needed to constrain the velocities, as follows. When using the velocity-less equations 
of motion, update the quaternions now using Equation (7.9) and the constraint is fully 
implemented. 

The velocity and angular rate are constrained in much the same way. 
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The impulse equation is used once again to distribute the velocity error between the two 
body’s velocities and angular rates. 
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The last step is to update the quaternions using the corrections computed in Equation (7.4). 


q“(r) = 6q‘'q“(r) 

q ^(0 = 6 q ^ q ^(0 
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As one may notice, the implementation of the pin constraint for the velocity-based equations of 
motion is about twice as expensive as the pin constraint in the velocity-less equations of 
motion. The pin constraint is also much more expensive than the very simple distance 
constraint. It is therefore important to try and use only particles and distance constraints 
whenever possible. Still, a single pin constraint connecting two rigid bodies involves fewer 
operations than assembling an entire network of eight particles and twelve distance constraints 
that do the same thing. 


8 The Angle Constraint 

Another very important constraint to implement is the angle constraint. It limits the amount one 
rigid body can rotate relative to another as shown in Figure 4. 



Figure 4 - Angle Constraint 

Since the equations of motion are written in terms of quaternions, and orientation limits are 
intuitively expressed as angular rotations around an axis, 0, a simple conversion is made. 
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The orientation of the second body relative to the first is computed in the first body’s coordinate 
frame: 


6q = q*(r)q“(r) 

SqL,=r(0Sqq‘'(0 


(8.3) 


The rotation into the first rigid body coordinate frame is necessary because quaternions are 
always expressed in the world coordinate frame. However, the angle limits are expressed as 
rotations relative to the first rigid body. It is now a simple matter to limit the rotations: 
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If none of the rotation limits have been exceeded, no further action is necessary. If an angle 
limit is violated, an orientation correction term is calculated in the world coordinate frame. 

= (8.5) 
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This correction quaternion, ([correction, must now be divided between the two rigid bodies in 
proportion to their inertias about the axis of the correction. 
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q“(r) = Slerp(q,q^„^^^^,„„,/“ )q“(r) 
q* (r) = Slerp (q, , /* ) q* (r) 

This involves quite a bit of work. To save time at the cost of accuracy, the magnitude of the 
inverse of the inertias can be used to compute f and and applied to Equation (8.9) instead. 
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When an angle limit is exceeded around any particular axis, the angular rates of the two rigid 
bodies around that axis must also be constrained. 
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A coefficient of restitution, s, is used to determine how fast the two bodies will rotate in the 
opposite direction after the limit is hit or whether they will simply stop rotating. The angular rate 
correction about any particular axis is set to zero if the orientation was not constrained around 
that axis. What is left of the angular rate correction is divided between the two rigid bodies 
using the same inertia ratios computed in Equations (8.8) or (8.10). 
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9 Collision Constraint 

Since the constraint iteration loop runs after the simulation loop (Section 10) it is possible for 
the effects of the constraints to overpower motion induced by some of the force models. This 
was observed in some of the early rag-doll simulations when parts of the rag-doll would 
penetrate into solid surfaces regardless of the collision model spring stiffness settings. The 
collision constraint solved this problem as the rag-doll’s collision constraints are iterated along 
with the rest of its pin constraints. The rag-doll application was the only one that necessitated 
the use of collision constraints instead of collision impulses and forces. 

As before, the collision detection algorithm returns a set of contact points, c„ surface normals, 
ft/, and the two rigid bodies involved in the collision. For each point and normal, the following 
algorithm is executed (all of the equations presented below should already be familiar). 
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At this point, the desired velocity after the collision is computed. 
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6x from Equation (9.1) is partitioned between the two rigid bodies using the impulse equation. 
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Similarly, 6v from Equation (9.6) is distributed between the two rigid bodies. 
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As before, the quaternions are updated last. 
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10 Simulation Loop 

The simulation loop is the main function that executes all of the force models, equations of 
motion and constraint solvers. The simulation runs at a minimum frequency for stability, so the 
entire simulation loop must be executed multiple times for larger step sizes. This minimum 
frequency depends on the force models used and is not a function of the constraints. In 
Q3.vr\es,f simulation IS typIcally 30 or 60 Hz. The number of iterations is computed using 
Equation (10.1). 
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As the number of constraint iterations is increased, the net residual error decreases and the 
constraints appear stronger. The number of iterations does not have to remain constant. It can 
be adjusted dynamically as a function of distance to the camera, for example. Conversely, the 
number of iterations can be decreased to make objects appear more flexible. For most 
applications, two or three iterations are sufficient. 

The ordering of the constraints is also critical. Every constraint that gets solved is affected by 
the constraints that follow, while the last constraint remains unaffected. To strengthen 
constraints, move them down in the list, to weaken constraints, move them up. Likewise, 
motion introduced by a force model can be overridden by a constraint. The simulation loop is 
illustrated in Figure 5. 
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Figure 5 - Simulation Loop 


11 Rag-doll Simulation 

The rag-doll simulation used in the game MX Unleashed, published by THQ, uses 18 rigid 
bodies, 17 pin, 17 angle and 9 collision constraints. The number of body parts can vary, but is 
usually determined by the needs of the animator. Spheres centered on some, but not all, of the 
body parts, were used to detect collisions. Depending on the performance needs of the 
application, spheres or the body parts themselves can be used. The setup used in MX 
Unleashed is illustrated in Figure 6. 









The positions of the pin constraints relative to the rigid bodies they connected were extracted 
from the geometry generated by the animator. The center of gravity was placed in the 
geometric center of each body part. The mass, inertias and angle limits are presented in 
Table 1. These mass properties are clearly not indicative of a real human body. Instead, they 
are settings that generated the desired motion in our game. The rag-doll simulation runs at 
30 Hz (one pass through the simulation) using two iterations of the constraint equations. 
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(135, 

0, 

0) 

Left Foot 

0.15 

(1, 1, 1) 

(-20, 

-20, 

0) 

(45, 

20, 

0) 

Right Upper Leg 

0.25 

(1, 1, 1) 

(-60, 

0, 

0) 

(0, 

0, 

45) 

Right Lower Leg 

0.20 

(1, 1, 1) 

(0, 

0, 

0) 

(135, 

0, 

0) 

Right Foot 

0.15 

(1, 1, 1) 

(-20, 

-20, 

0) 

(45, 

20, 

0) 


Table 1 - Rag-doll Mass Properties and Angle Constraints 












































12 Performance 


The rag-doll simulation described in Section 1 1 was run on both the PlayStation 2 and XBox. 
On the PS2, the equations of motion and the constraint equations run entirely on the vector 
unit (VUO) in macromode. The simulation loop runs on the CPU. On the XBox, the quaternion 
math was hand optimized using the SSE (Streaming SIMD Extensions) instruction set but the 
compiler optimized the rest of the code. The timing results are listed in Table 2. The collision 
detection times are not included in these results. However, the collision constraint times 
include a small percentage of time that should actually be attributed to the collision detection 
algorithm. 


Algorithm 

Number of 
Calls 

Total Time 
PS2 
(mSec) 

Time Per Call 
PS2 
(mSec) 

Total Time 
XBox 

(mSec) 

Time Per Call 
XBox 
(mSec) 

Pin Constraint 

34 

(2x 17) 

0.17 

0.0051 

0.15 

0.0044 

Angle Constraint 

34 

(2x 17) 

0.15 

0.0045 

0.21 

0.0062 

Collision Constraint 

18 

(2x9) 

0.15 

0.0086 

0.13 

0.0072 

Equations of Motion 

18 

(1 X 18) 

0.03 

0.0017 

0.06 

0.0033 

Total 


0.50 


0.55 



Table 2 - Rag-doll Performance on the PlayStation 2 and XBox 


13 Conclusion 

The Verlet-like integrator made it possible to integrate systems of sufficient stiffness for our 
applications. This removed the need to use implicit integrators which can be difficult to 
program and expensive to solve. Penalty methods could then be used to resolve resting forces 
instead of having to solve for them explicitly. The difficulty in getting rigid bodies to sit still on a 
surface when using the impulse method was also avoided. 

The addition of the rotational degrees of freedom simplified the simulation of rigid bodies and 
made more traditional pin and angle constraints possible. Replacing the last position and 
quaternion with a velocity and angular rate solved the problem of computing damping forces 
based on a noisy difference equation. However, this made the constraint equations more 
expensive as a result. 

By trading the stability associated with implicit integration for performance, we feel that this 
simulation is the ideal trade-off for this generation of console games. As memory and CPU 
speeds increase, it may not be as critical to keep the code as fast as possible and implicit 
integration using large sparse matrices may become the best solution for the future. 


Please feel free to contact the authors should you have any questions, corrections, 
suggestions or improvements. 
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Appendix A - Quaternion Operations 

Basic quaternion: 



Unit quaternion; 





Multiplication; 



First Derivative: 


q(q,e> =-q 

2 CO. 


Second Derivative: 


-2q q 


■ ■ j • \ 1 

q(q,e>,e>) = -q 

2 CO., 


Conversion to Angle Axis Vector: 


ComputeAngleAxis(q) = 


where 


0 = 2cos ' 


Conversion from Angle Axis Vector: 


A 0 

0 cos — 
2 

. 0 

, v^sin — 

1 ^ 2 

ComputeQuaternion(v) = — 

0 . 0 
v„ sin — 

^ 2 

. 0 
sin — 

L 2_ 


where 


0 =|v| 



Spherical Linear Interpolation: 


where 


Slerp(q“,q\r) 


aql + P^* 


sin((l-r)0) 

a = — - 

sin0 

sii^ 

sin0 

0 = cos ' (q" -q*^ 



